{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "append-SMB-ATL11\n",
    "================\n",
    "\n",
    "Interpolates modeled firn estimates to the coordinates of an ATL11 file\n",
    "\n",
    "#### Python Dependencies\n",
    "- [numpy: Scientific Computing Tools For Python](https://numpy.org)  \n",
    "- [netCDF4: Python interface to the netCDF C library](https://unidata.github.io/netcdf4-python/netCDF4/index.html)  \n",
    "- [h5py: Python interface for Hierarchal Data Format 5 (HDF5)](https://www.h5py.org/)  \n",
    "- [pointCollection: Utilities for organizing and manipulating point data](https://github.com/SmithB/pointCollection)  \n",
    "\n",
    "This notebook uses Jupyter widgets to set parameters for calculating the firn estimates.  \n",
    "The widgets can be installed as described below.  \n",
    "```\n",
    "pip3 install --user ipywidgets\n",
    "jupyter nbextension install --user --py widgetsnbextension\n",
    "jupyter nbextension enable --user --py widgetsnbextension\n",
    "jupyter-notebook\n",
    "```\n",
    "\n",
    "#### Load modules "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import re\n",
    "import SMBcorr\n",
    "import h5py\n",
    "import numpy as np\n",
    "import pointCollection as pc\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Choose parameters to run\n",
    "- Input ATL11 file  \n",
    "- Directory containing SMB models  \n",
    "- SMB region to run  \n",
    "- SMB model to run  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose file to read and append\n",
    "value=''\n",
    "inputText = widgets.Text(\n",
    "    value=value,\n",
    "    description='File:',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "# text menu to set model base directory\n",
    "value=os.getcwd()\n",
    "directoryText = widgets.Text(\n",
    "    value=value,\n",
    "    description='Directory:',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "# dropdown menu for setting region\n",
    "regionDropdown = widgets.Dropdown(\n",
    "    options=['AA','GL'],\n",
    "    value='GL',\n",
    "    description='Region:',\n",
    "    disabled=False,\n",
    ")\n",
    "\n",
    "# dropdown menu for setting SMB model\n",
    "modelDropdown = widgets.Dropdown(\n",
    "    options=['MAR','RACMO','MERRA2-hybrid'],\n",
    "    value='MAR',\n",
    "    description='Model:',\n",
    "    disabled=False,\n",
    ")\n",
    "\n",
    "# dropdown menu for setting model version\n",
    "models = dict(AA={}, GL={})\n",
    "# MAR\n",
    "models['GL']['MAR'] = []\n",
    "# models['GL']['MAR'].append('MARv3.9-ERA')\n",
    "# models['GL']['MAR'].append('MARv3.10-ERA')\n",
    "# models['GL']['MAR'].append('MARv3.11-NCEP')\n",
    "models['GL']['MAR'].append('MARv3.11-ERA')\n",
    "models['GL']['MAR'].append('MARv3.11.2-ERA-7.5km')\n",
    "models['GL']['MAR'].append('MARv3.11.2-ERA-10km')\n",
    "models['GL']['MAR'].append('MARv3.11.2-ERA-15km')\n",
    "models['GL']['MAR'].append('MARv3.11.2-ERA-20km')\n",
    "models['GL']['MAR'].append('MARv3.11.2-NCEP-20km')\n",
    "# RACMO\n",
    "models['GL']['RACMO'] = []\n",
    "# models['GL']['RACMO'].append('RACMO2.3-XGRN11')\n",
    "# models['GL']['RACMO'].append('RACMO2.3p2-XGRN11')\n",
    "models['GL']['RACMO'].append('RACMO2.3p2-FGRN055')\n",
    "# MERRA2-hybrid\n",
    "models['GL']['MERRA2-hybrid'] = []\n",
    "# models['GL']['MERRA2-hybrid'].append('GSFC-fdm-v0')\n",
    "models['GL']['MERRA2-hybrid'].append('GSFC-fdm-v1')\n",
    "models['AA']['MERRA2-hybrid'] = []\n",
    "# models['AA']['MERRA2-hybrid'].append('GSFC-fdm-v0')\n",
    "models['AA']['MERRA2-hybrid'].append('GSFC-fdm-v1')\n",
    "versionDropdown = widgets.Dropdown(\n",
    "    options=models[regionDropdown.value][modelDropdown.value],\n",
    "    value=models[regionDropdown.value][modelDropdown.value][0],\n",
    "    description='Version:',\n",
    "    disabled=False,\n",
    ")\n",
    "\n",
    "# function for updating the model version\n",
    "def set_version(sender):\n",
    "    versionDropdown.options=models[regionDropdown.value][modelDropdown.value]\n",
    "    versionDropdown.value=models[regionDropdown.value][modelDropdown.value][0]\n",
    "\n",
    "# watch widgets for changes\n",
    "regionDropdown.observe(set_version)\n",
    "modelDropdown.observe(set_version)\n",
    "\n",
    "# display widgets for setting parameters\n",
    "widgets.VBox([inputText,directoryText,regionDropdown,\n",
    "    modelDropdown,versionDropdown])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function to convert ATL11 time variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PURPOSE: convert time from delta seconds into Julian and year-decimal\n",
    "def convert_delta_time(delta_time, gps_epoch=1198800018.0):\n",
    "    # calculate gps time from delta_time\n",
    "    gps_seconds = gps_epoch + delta_time\n",
    "    time_leaps = SMBcorr.time.count_leap_seconds(gps_seconds)\n",
    "    # calculate julian time\n",
    "    time_julian = 2400000.5 + SMBcorr.time.convert_delta_time(gps_seconds - time_leaps,\n",
    "        epoch1=(1980,1,6,0,0,0), epoch2=(1858,11,17,0,0,0), scale=1.0/86400.0)\n",
    "    # convert to calendar date with convert_julian.py\n",
    "    Y,M,D,h,m,s = SMBcorr.convert_julian(time_julian,FORMAT='tuple')\n",
    "    # calculate year-decimal time\n",
    "    time_decimal = SMBcorr.convert_calendar_decimal(Y,M,DAY=D,HOUR=h,MINUTE=m,SECOND=s)\n",
    "    # return both the Julian and year-decimal formatted dates\n",
    "    return dict(julian=time_julian, decimal=time_decimal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function to set the coordinate projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PURPOSE: set the projection parameters based on the region name\n",
    "def set_projection(REGION):\n",
    "    if (REGION == 'AA'):\n",
    "        projection_flag = 'EPSG:3031'\n",
    "    elif (REGION == 'GL'):\n",
    "        projection_flag = 'EPSG:3413'\n",
    "    return projection_flag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run interpolation program for model and ATL11 file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read input file\n",
    "input_file = os.path.expanduser(inputText.value)\n",
    "field_dict = {None:('delta_time','h_corr','x','y')}\n",
    "D11 = pc.data().from_h5(input_file, field_dict=field_dict)\n",
    "# check if running crossover or along-track ATL11\n",
    "if (D11.h_corr.ndim == 3):\n",
    "    nseg,ncycle,ncross = D11.shape\n",
    "else:\n",
    "    nseg,ncycle = D11.shape\n",
    "        \n",
    "# get projection of input coordinates\n",
    "EPSG = set_projection(regionDropdown.value)\n",
    "\n",
    "# extract parameters from widgets\n",
    "base_dir=os.path.expanduser(directoryText.value)\n",
    "if (modelDropdown.value == 'MAR'):\n",
    "    match_object=re.match('(MARv\\d+\\.\\d+(.\\d+)?)',versionDropdown.value)\n",
    "    MAR_VERSION=match_object.group(0)\n",
    "    MAR_REGION=dict(GL='Greenland',AA='Antarctic')[regionDropdown.value]\n",
    "    # model subdirectories\n",
    "    SUBDIRECTORY=dict(AA={}, GL={})\n",
    "    SUBDIRECTORY['GL']['MARv3.9-ERA']=['ERA_1958-2018_10km','daily_10km']\n",
    "    SUBDIRECTORY['GL']['MARv3.10-ERA']=['ERA_1958-2019-15km','daily_15km']\n",
    "    SUBDIRECTORY['GL']['MARv3.11-NCEP']=['NCEP1_1948-2020_20km','daily_20km']\n",
    "    SUBDIRECTORY['GL']['MARv3.11-ERA']=['ERA_1958-2019-15km','daily_15km']\n",
    "    SUBDIRECTORY['GL']['MARv3.11.2-ERA-7.5km']=['7.5km_ERA5']\n",
    "    SUBDIRECTORY['GL']['MARv3.11.2-ERA-10km']=['10km_ERA5']\n",
    "    SUBDIRECTORY['GL']['MARv3.11.2-ERA-15km']=['15km_ERA5']\n",
    "    SUBDIRECTORY['GL']['MARv3.11.2-ERA-20km']=['20km_ERA5']\n",
    "    SUBDIRECTORY['GL']['MARv3.11.2-NCEP-20km']=['20km_NCEP1']\n",
    "    MAR_MODEL=SUBDIRECTORY[regionDropdown.value][versionDropdown.value]\n",
    "    DIRECTORY=os.path.join(base_dir,'MAR',MAR_VERSION,MAR_REGION,*MAR_MODEL)    \n",
    "    # variable coordinates\n",
    "    KWARGS=dict(AA={}, GL={})\n",
    "    KWARGS['GL']['MARv3.9-ERA'] = dict(XNAME='X10_153',YNAME='Y21_288')\n",
    "    KWARGS['GL']['MARv3.10-ERA'] = dict(XNAME='X10_105',YNAME='Y21_199')\n",
    "    KWARGS['GL']['MARv3.11-NCEP'] = dict(XNAME='X12_84',YNAME='Y21_155')\n",
    "    KWARGS['GL']['MARv3.11-ERA'] = dict(XNAME='X10_105',YNAME='Y21_199')\n",
    "    KWARGS['GL']['MARv3.11.2-ERA-7.5km'] = dict(XNAME='X12_203',YNAME='Y20_377')\n",
    "    KWARGS['GL']['MARv3.11.2-ERA-10km'] = dict(XNAME='X10_153',YNAME='Y21_288')\n",
    "    KWARGS['GL']['MARv3.11.2-ERA-15km'] = dict(XNAME='X10_105',YNAME='Y21_199')\n",
    "    KWARGS['GL']['MARv3.11.2-ERA-20km'] = dict(XNAME='X12_84',YNAME='Y21_155')\n",
    "    KWARGS['GL']['MARv3.11.2-NCEP-20km'] = dict(XNAME='X12_84',YNAME='Y21_155')\n",
    "    MAR_KWARGS=KWARGS[regionDropdown.value][versionDropdown.value]\n",
    "    # output variable keys for both direct and derived fields\n",
    "    KEYS = ['zsurf','zfirn','zmelt','zsmb','zaccum']\n",
    "    # HDF5 longname attributes for each variable\n",
    "    LONGNAME = {}\n",
    "    LONGNAME['zsurf'] = \"Snow Height Change\"\n",
    "    LONGNAME['zfirn'] = \"Snow Height Change due to Compaction\"\n",
    "    LONGNAME['zmelt'] = \"Snow Height Change due to Surface Melt\"\n",
    "    LONGNAME['zsmb'] = \"Snow Height Change due to Surface Mass Balance\"\n",
    "    LONGNAME['zaccum'] = \"Snow Height Change due to Surface Accumulation\"\n",
    "elif (modelDropdown.value == 'RACMO'):\n",
    "    RACMO_VERSION,RACMO_MODEL=versionDropdown.value.split('-')\n",
    "    # output variable keys\n",
    "    KEYS = ['zsurf']\n",
    "    # HDF5 longname attributes for each variable\n",
    "    LONGNAME = {}\n",
    "    LONGNAME['zsurf'] = \"Snow Height Change\"\n",
    "elif (modelDropdown.value == 'MERRA2-hybrid'):\n",
    "    MERRA2_VERSION,=re.findall('GSFC-fdm-(.*?)$',versionDropdown.value)\n",
    "    # MERRA-2 hybrid directory\n",
    "    DIRECTORY=os.path.join(base_dir,'MERRA2_hybrid',MERRA2_VERSION)\n",
    "    MERRA2_REGION = dict(AA='ais',GL='gris')[regionDropdown.value]\n",
    "    # output variable keys for both direct and derived fields\n",
    "    KEYS = ['zsurf','zfirn','zsmb']\n",
    "    # HDF5 longname attributes for each variable\n",
    "    LONGNAME = {}\n",
    "    LONGNAME['zsurf'] = \"Snow Height Change\"\n",
    "    LONGNAME['zfirn'] = \"Snow Height Change due to Compaction\"\n",
    "    LONGNAME['zsmb'] = \"Snow Height Change due to Surface Mass Balance\"\n",
    "    \n",
    "# check if running crossover or along track\n",
    "if (D11.h_corr.ndim == 3):\n",
    "    # allocate for output height for crossover data \n",
    "    OUTPUT = {}\n",
    "    for key in KEYS:\n",
    "        OUTPUT[key] = np.ma.zeros((nseg,ncycle,ncross),fill_value=np.nan)\n",
    "        OUTPUT[key].mask = np.ones((nseg,ncycle,ncross),dtype=np.bool)\n",
    "        OUTPUT[key].interpolation = np.zeros((nseg,ncycle,ncross),dtype=np.uint8)\n",
    "    # for each cycle of ICESat-2 ATL11 data\n",
    "    for c in range(ncycle):\n",
    "        # check that there are valid crossovers\n",
    "        cross = [xo for xo in range(ncross) if\n",
    "            np.any(np.isfinite(D11.delta_time[:,c,xo]))]\n",
    "        # for each valid crossing\n",
    "        for xo in cross:\n",
    "            # find valid crossovers\n",
    "            i, = np.nonzero(np.isfinite(D11.delta_time[:,c,xo]))\n",
    "            # convert from delta time to decimal-years\n",
    "            tdec = convert_delta_time(D11.delta_time[i,c,xo])['decimal']\n",
    "            if (modelDropdown.value == 'MAR'):\n",
    "                # read and interpolate daily MAR outputs\n",
    "                ZN4 = SMBcorr.interpolate_mar_daily(DIRECTORY, EPSG,\n",
    "                    MAR_VERSION, tdec, D11.x[i,c,xo], D11.y[i,c,xo],\n",
    "                    VARIABLE='ZN4', SIGMA=1.5, FILL_VALUE=np.nan, **MAR_KWARGS)\n",
    "                ZN5 = SMBcorr.interpolate_mar_daily(DIRECTORY, EPSG,\n",
    "                    MAR_VERSION, tdec, D11.x[i,c,xo], D11.y[i,c,xo],\n",
    "                    VARIABLE='ZN5', SIGMA=1.5, FILL_VALUE=np.nan, **MAR_KWARGS)\n",
    "                ZN6 = SMBcorr.interpolate_mar_daily(DIRECTORY, EPSG,\n",
    "                    MAR_VERSION, tdec, D11.x[i,c,xo], D11.y[i,c,xo],\n",
    "                    VARIABLE='ZN6', SIGMA=1.5, FILL_VALUE=np.nan, **MAR_KWARGS)\n",
    "                # set attributes to output for iteration\n",
    "                OUTPUT['zfirn'].data[i,c,xo] = np.copy(ZN4.data)\n",
    "                OUTPUT['zfirn'].mask[i,c,xo] = np.copy(ZN4.mask)\n",
    "                OUTPUT['zfirn'].interpolation[i,c,xo] = np.copy(ZN4.interpolation)\n",
    "                OUTPUT['zsurf'].data[i,c,xo] = np.copy(ZN6.data)\n",
    "                OUTPUT['zsurf'].mask[i,c,xo] = np.copy(ZN6.mask)\n",
    "                OUTPUT['zsurf'].interpolation[i,c,xo] = np.copy(ZN6.interpolation)\n",
    "                OUTPUT['zmelt'].data[i,c,xo] = np.copy(ZN5.data)\n",
    "                OUTPUT['zmelt'].mask[i,c,xo] = np.copy(ZN5.mask)\n",
    "                OUTPUT['zmelt'].interpolation[i,c,xo] = np.copy(ZN5.interpolation)\n",
    "                # calculate derived fields\n",
    "                OUTPUT['zsmb'].data[i,c,xo] = ZN6.data - ZN4.data\n",
    "                OUTPUT['zsmb'].mask[i,c,xo] = ZN4.mask | ZN6.mask\n",
    "                OUTPUT['zaccum'].data[i,c,xo] = ZN6.data - ZN4.data - ZN5.data\n",
    "                OUTPUT['zaccum'].mask[i,c,xo] = ZN4.mask | ZN5.mask | ZN6.mask\n",
    "            elif (modelDropdown.value == 'RACMO'):\n",
    "                # read and interpolate daily RACMO outputs\n",
    "                hgtsrf = SMBcorr.interpolate_racmo_daily(base_dir, EPSG,\n",
    "                    RACMO_MODEL, tdec, D11.x[i,c,xo], D11.y[i,c,xo],\n",
    "                    VARIABLE='hgtsrf', SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "                # set attributes to output for iteration\n",
    "                OUTPUT['zsurf'].data[i,c,xo] = np.copy(hgtsrf.data)\n",
    "                OUTPUT['zsurf'].mask[i,c,xo] = np.copy(hgtsrf.mask)\n",
    "                OUTPUT['zsurf'].interpolation[i,c,xo] = np.copy(hgtsrf.interpolation)\n",
    "            elif (modelDropdown.value == 'MERRA2-hybrid'):\n",
    "                # read and interpolate 5-day MERRA2-Hybrid outputs\n",
    "                FAC = SMBcorr.interpolate_merra_hybrid(DIRECTORY, EPSG,\n",
    "                    MERRA2_REGION, tdec, D11.x[i,c,xo], D11.y[i,c,xo],\n",
    "                    VERSION=MERRA2_VERSION, VARIABLE='FAC',\n",
    "                    SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "                smb = SMBcorr.interpolate_merra_hybrid(DIRECTORY, EPSG,\n",
    "                    MERRA2_REGION, tdec, D11.x[i,c,xo], D11.y[i,c,xo],\n",
    "                    VERSION=MERRA2_VERSION, VARIABLE='cum_smb_anomaly',\n",
    "                    SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "                height = SMBcorr.interpolate_merra_hybrid(DIRECTORY, EPSG,\n",
    "                    MERRA2_REGION, tdec, D11.x[i,c,xo], D11.y[i,c,xo],\n",
    "                    VERSION=MERRA2_VERSION, VARIABLE='height',\n",
    "                    SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "                # set attributes to output for iteration\n",
    "                OUTPUT['zfirn'].data[i,c,xo] = np.copy(FAC.data)\n",
    "                OUTPUT['zfirn'].mask[i,c,xo] = np.copy(FAC.mask)\n",
    "                OUTPUT['zfirn'].interpolation[i,c,xo] = np.copy(FAC.interpolation)\n",
    "                OUTPUT['zsurf'].data[i,c,xo] = np.copy(height.data)\n",
    "                OUTPUT['zsurf'].mask[i,c,xo] = np.copy(height.mask)\n",
    "                OUTPUT['zsurf'].interpolation[i,c,xo] = np.copy(height.interpolation)\n",
    "                OUTPUT['zsmb'].data[i,c,xo] = np.copy(smb.data)\n",
    "                OUTPUT['zsmb'].mask[i,c,xo] = np.copy(smb.mask)\n",
    "                OUTPUT['zsmb'].interpolation[i,c,xo] = np.copy(smb.interpolation)\n",
    "else:        \n",
    "    # allocate for output height for along-track data \n",
    "    OUTPUT = {}\n",
    "    for key in KEYS:\n",
    "        OUTPUT[key] = np.ma.zeros((nseg,ncycle),fill_value=np.nan)\n",
    "        OUTPUT[key].mask = np.ones((nseg,ncycle),dtype=np.bool)\n",
    "        OUTPUT[key].interpolation = np.zeros((nseg,ncycle),dtype=np.uint8)    \n",
    "    # check that there are valid elevations\n",
    "    cycle = [c for c in range(ncycle) if\n",
    "        np.any(np.isfinite(D11.delta_time[:,c]))]\n",
    "    # for each valid cycle of ICESat-2 ATL11 data\n",
    "    for c in cycle:\n",
    "        # find valid elevations\n",
    "        i, = np.nonzero(np.isfinite(D11.delta_time[:,c]))\n",
    "        # convert from delta time to decimal-years\n",
    "        tdec = convert_delta_time(D11.delta_time[i,c])['decimal']\n",
    "        if (modelDropdown.value == 'MAR'):\n",
    "            # read and interpolate daily MAR outputs\n",
    "            ZN4 = SMBcorr.interpolate_mar_daily(DIRECTORY, EPSG,\n",
    "                MAR_VERSION, tdec, D11.x[i,c], D11.y[i,c], VARIABLE='ZN4',\n",
    "                SIGMA=1.5, FILL_VALUE=np.nan, **MAR_KWARGS)\n",
    "            ZN5 = SMBcorr.interpolate_mar_daily(DIRECTORY, EPSG,\n",
    "                MAR_VERSION, tdec, D11.x[i,c], D11.y[i,c], VARIABLE='ZN5',\n",
    "                SIGMA=1.5, FILL_VALUE=np.nan, **MAR_KWARGS)\n",
    "            ZN6 = SMBcorr.interpolate_mar_daily(DIRECTORY, EPSG,\n",
    "                MAR_VERSION, tdec, D11.x[i,c], D11.y[i,c], VARIABLE='ZN6',\n",
    "                SIGMA=1.5, FILL_VALUE=np.nan, **MAR_KWARGS)\n",
    "            # set attributes to output for iteration\n",
    "            OUTPUT['zfirn'].data[i,c] = np.copy(ZN4.data)\n",
    "            OUTPUT['zfirn'].mask[i,c] = np.copy(ZN4.mask)\n",
    "            OUTPUT['zfirn'].interpolation[i,c] = np.copy(ZN4.interpolation)\n",
    "            OUTPUT['zsurf'].data[i,c] = np.copy(ZN6.data)\n",
    "            OUTPUT['zsurf'].mask[i,c] = np.copy(ZN6.mask)\n",
    "            OUTPUT['zsurf'].interpolation[i,c] = np.copy(ZN6.interpolation)\n",
    "            OUTPUT['zmelt'].data[i,c] = np.copy(ZN5.data)\n",
    "            OUTPUT['zmelt'].mask[i,c] = np.copy(ZN5.mask)\n",
    "            OUTPUT['zmelt'].interpolation[i,c] = np.copy(ZN5.interpolation)\n",
    "            # calculate derived fields\n",
    "            OUTPUT['zsmb'].data[i,c] = ZN6.data - ZN4.data\n",
    "            OUTPUT['zsmb'].mask[i,c] = ZN4.mask | ZN6.mask\n",
    "            OUTPUT['zaccum'].data[i,c] = ZN6.data - ZN4.data - ZN5.data\n",
    "            OUTPUT['zaccum'].mask[i,c] = ZN4.mask | ZN5.mask | ZN6.mask            \n",
    "        elif (modelDropdown.value == 'RACMO'):\n",
    "            # read and interpolate daily RACMO outputs\n",
    "            hgtsrf = SMBcorr.interpolate_racmo_daily(base_dir, EPSG,\n",
    "                RACMO_MODEL, tdec, D11.x[i,c], D11.y[i,c],\n",
    "                VARIABLE='hgtsrf', SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "            # set attributes to output for iteration\n",
    "            OUTPUT['zsurf'].data[i,c] = np.copy(hgtsrf.data)\n",
    "            OUTPUT['zsurf'].mask[i,c] = np.copy(hgtsrf.mask)\n",
    "            OUTPUT['zsurf'].interpolation[i,c] = np.copy(hgtsrf.interpolation)\n",
    "        elif (modelDropdown.value == 'MERRA2-hybrid'):\n",
    "            # read and interpolate 5-day MERRA2-Hybrid outputs\n",
    "            FAC = SMBcorr.interpolate_merra_hybrid(DIRECTORY, EPSG,\n",
    "                MERRA2_REGION, tdec, D11.x[i,c], D11.y[i,c],\n",
    "                VERSION=MERRA2_VERSION, VARIABLE='FAC',\n",
    "                SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "            smb = SMBcorr.interpolate_merra_hybrid(DIRECTORY, EPSG,\n",
    "                MERRA2_REGION, tdec, D11.x[i,c], D11.y[i,c],\n",
    "                VERSION=MERRA2_VERSION, VARIABLE='cum_smb_anomaly',\n",
    "                SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "            height = SMBcorr.interpolate_merra_hybrid(DIRECTORY, EPSG,\n",
    "                MERRA2_REGION, tdec, D11.x[i,c], D11.y[i,c],\n",
    "                VERSION=MERRA2_VERSION, VARIABLE='height',\n",
    "                SIGMA=1.5, FILL_VALUE=np.nan)\n",
    "            # set attributes to output for iteration\n",
    "            OUTPUT['zfirn'].data[i,c] = np.copy(FAC.data)\n",
    "            OUTPUT['zfirn'].mask[i,c] = np.copy(FAC.mask)\n",
    "            OUTPUT['zfirn'].interpolation[i,c] = np.copy(FAC.interpolation)\n",
    "            OUTPUT['zsurf'].data[i,c] = np.copy(height.data)\n",
    "            OUTPUT['zsurf'].mask[i,c] = np.copy(height.mask)\n",
    "            OUTPUT['zsurf'].interpolation[i,c] = np.copy(height.interpolation)\n",
    "            OUTPUT['zsmb'].data[i,c] = np.copy(smb.data)\n",
    "            OUTPUT['zsmb'].mask[i,c] = np.copy(smb.mask)\n",
    "            OUTPUT['zsmb'].interpolation[i,c] = np.copy(smb.interpolation)\n",
    "                \n",
    "# verify mask values\n",
    "for key in KEYS:\n",
    "    OUTPUT[key].mask |= (OUTPUT[key].data == OUTPUT[key].fill_value) | \\\n",
    "        np.isnan(OUTPUT[key].data)\n",
    "    OUTPUT[key].data[OUTPUT[key].mask] = OUTPUT[key].fill_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create plot showing elevation change and total modeled height change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (D11.h_corr.ndim == 2):\n",
    "    fig,(ax1,ax2) = plt.subplots(ncols=2)\n",
    "    c1,c2 = (0,1)\n",
    "    i, = np.nonzero(np.isfinite(D11.h_corr[:,c1]) & np.isfinite(D11.h_corr[:,c2]) & \n",
    "        (~OUTPUT['zsurf'].mask[:,c1]) & (~OUTPUT['zsurf'].mask[:,c2]))\n",
    "    # height change from ICESat-2\n",
    "    ax1.scatter(D11.x[i,c1],D11.y[i,c1],\n",
    "        c=D11.h_corr[i,c2]-D11.h_corr[i,c1],\n",
    "        vmin=-5,vmax=0.5,s=0.5)\n",
    "    ax1.set_title = 'ATL11'\n",
    "    # firn height change\n",
    "    ax2.scatter(D11.x[i,c1],D11.y[i,c1],\n",
    "        c=OUTPUT['zsurf'][i,c2]-OUTPUT['zsurf'][i,c1],\n",
    "        vmin=-5,vmax=0.5,s=0.5)\n",
    "    ax2.set_title = versionDropdown.value\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Append data to input ATL11 HDF5 file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# append input HDF5 file with new firn model outputs\n",
    "fileID = h5py.File(os.path.expanduser(inputText.value),'a')\n",
    "fileID.create_group(versionDropdown.value)\n",
    "h5 = {}\n",
    "for key in KEYS:\n",
    "    val = '{0}/{1}'.format(versionDropdown.value,key)\n",
    "    h5[key] = fileID.create_dataset(val, OUTPUT[key].shape,\n",
    "        data=OUTPUT[key], dtype=OUTPUT[key].dtype,\n",
    "        compression='gzip', fillvalue=OUTPUT[key].fill_value)\n",
    "    h5[key].attrs['units'] = \"m\"\n",
    "    h5[key].attrs['long_name'] = LONGNAME[key]\n",
    "    h5[key].attrs['coordinates'] = \"../delta_time ../latitude ../longitude\"\n",
    "    h5[key].attrs['model'] = versionDropdown.value\n",
    "# close the output HDF5 file\n",
    "fileID.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
